A Universal Map for Fractal Structures in Weak Solitary Wave Interactions 
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Fractal scatterings in weak solitary wave interactions is analyzed for generalized nonlinear 
Schrodiger equations (GNLS). Using asymptotic methods, these weak interactions are reduced to 
a universal second-order map. This map gives the same fractal scattering patterns as those in the 
GNLS equations both qualitatively and quantitatively. Scaling laws of these fractals are also derived. 
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PACS numbers: 42.65.Tg, 05.45.Yv, 42.81. Dp 

Solitary wave interactions is a fascinating mathemati- 
cal phenomenon, and it arises in numerous physical appli- 
cations such as water waves and nonlinear optics ^1, 2, 3]. 
Strong interactions occur when two solitary waves are ini- 
tially far apart but move toward each other at moderate 
or large speeds. Weak interactions would occur if the two 
waves are initially well separated, and their relative veloc- 
ities are small. In integrable wave equations, strong inter- 
actions of solitary waves are elastic, and their weak inter- 
actions exhibit interesting but simple dynamics [H, 0, 0] ■ 
In non-integrable systems, however, solitary wave inter- 
actions can be extremely complicated. Indeed, one of the 
most important developments in the nonlinear wave the- 
ory in recent years is the discovery of fractal scatterings 
of solitary wave interactions in non-integrable equations 
d, H, 0, B S [HI- On the analysis of fractal scat- 
terings, some progress has been made. For strong inter- 
actions, various collective-coordinate ODE models based 
on qualitative variational methods have been derived and 
analyzed [1, 0, H, [HI • From the variational ODEs, a sepa- 
ratrix map was derived, showing chaotic scatterings [11]. 
On weak interactions, a simple asymptotically-accurate 
ODE model was derived for the generalized NLS equa- 
tions This ODE system offered the first glimpse of 
universal fractal scatterings in weak wave interactions, 
but these fractal patterns were not analyzed. 

In this letter, we analyze fractal scattering patterns 
in weak solitary wave interactions for generalized nonlin- 
ear Schrodiger equations. Using asymptotic methods, we 
reduce these weak interactions to a simple second-order 
map which contains no free parameters. It is shown that 
this universal map gives a complete characterization of 
fractal structures in these wave interactions. In addition, 
the scaling laws of these fractals for different initial con- 
ditions are analytically derived. These results provide 
a deep understanding of weak solitary-wave interactions 
for various physical applications. 

The generalized nonlinear Schrodingcr equations we 
consider in this paper are 



iUt + U^., + F{\U\^)U = 0, 



(1) 



where F{-) is a general function. These equations gov- 
ern various physical wave phenomena in nonlinear a 
tics, fiber communications and fluid dynamics P, 



This equation admits solitary waves of the form U = 
^{x — O^*"^: where $(6') is a localized positive func- 
tion, ^ = Vt + xq is the wave's center position, and 
(t> = \V{x — + (/5 + 3^^)^ ~ ^0 is the wave's phase. 
This wave has four free parameters: velocity V , ampli- 
tude parameter (3, initial position xq, and initial phase 770. 
In weak interactions, two such solitary waves are initially 
well separated and have small relative velocities and am- 
plitude differences. Then they would interfere with each 
other through tail overlapping. When time goes to infin- 
ity, they cither separate from each other or form a bound 
state. The exit velocity, defined as AI/qq = 1^2 — ^i|t->ooj 
depends on the initial conditions of the two waves. 

To study weak interactions in Eq. ([T]), we select 
two different nonlinearities which are cubic-quintic and 
quadratic-cubic respectively: 

F{\U\^)^\U\' + ^\U\\ F(|C/|2) = |C/|2 + <5|[/|. (2) 

Here 7 and 5 are real parameters. To illustrate results of 
weak interactions, we take 7 = 0.0003 and 5 — —0.0015. 
The initial conditions are taken as 

= -xq,2 = -5, y^fc = 0, /?o,fc = 1 (fc = 1, 2), (3) 

'/'o.i = 0, and the initial phase difference A0o = 4'o,2 — 
(/>o,i is used as the control parameter. In our numerical 
simulations of Eq. ([T]), the discrete Fourier transform 
is used to evaluate the spatial derivative dxx, while the 
fourth-order Runge-Kutta method is used to advance in 
time. The exit velocity A14o versus A0o graphs for these 
two nonlinearities are plotted in Fig. [T] These graphs are 
fractals [l3| • It is amazing that these fractals appear for 
such small values of 7 and 5, where Eq. ([1]) is simply 
a weakly perturbed NLS equation. Notice that the frac- 
tals for these two different nonlinearities are very similar, 
signaling their universality in weak wave interactions. 

To analyze this fractal-scattering phenomena, the 
Karpman-Solov'ev method [3| was applied, and the fol- 
lowing simple set of dynamical equations for solitary wave 
parameters were derived |10]: 

cosV'e'', Vrr = (1 + e) sin^/ie^. 



Crr 



(4) 



Here = A0, C = -VP^L 



T = i/l6/33/2c2/F t, e = Pl{2pPp) - 1, 



(5) 
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FIG. 1: Exit velocity versus initial phase difference graphs 
for initial conditions ([31 : (a) cubic-quintic nonlinearity; (b) 
quadratic-cubic nonlinearity. 



and A0 are the distance and phase difference between 
the two waves, (3 = (/?i,o + /32,o)/2, c is the tail coefficient 
of the solitary wave with propagation constant /3, and 
P{f3) is the power function of the wave. These reduced 
ODEs capture fractal scatterings of weak wave interac- 
tions such as in Fig. 1 both qualitatively and quanti- 
tatively and they represent an important first step 
toward the understanding of these phenomena. However, 
the fractal structures in the PDEs IJ) and ODEs ([U have 
not been analyzed previously. Below we give a complete 
characterization for the first time of this fractal scatter- 
ing by analyzing the ODEs Q. 

If e = 0, Eq. (|4]) is integrable. It has two conserved 
quantities, energy E and momentum M: 

E^{(^ - V?)/2 - cos ?/>, Af = CV- - e« sin ^. (6) 

Introducing two complex quantities 



C = ^{E + tM)/2, F = -acoth[(C + iyj)/2C)]/C, (7) 

the analytical solution of Eq. ^ can be found to be 

Y{t) = In [2C^csch'CoiT - tq + Fq)] , (8) 

where Y = C + iip, and Cq, Fq are the initial values 
of C and F. The third conserved quantity of Eq. (j4]) 
is Im(i^). Behaviors of the above integrable solutions 
should be noted. When > 0, or £; < but M ^ 0, 
C — oo (a degenerate saddle point) as r ^ oo, and thus 
these solutions are escape orbits. When E < and M = 
0, the orbits are periodic with period Tp — / \/\E\. 
Orbits with E = M = Q separate the escape orbits from 
the periodic ones, hence we call them separatrix orbits. 
The formulae for separatrix orbits can be readily found 
to be 



(9) 



Here A = e'^^' , C-A/ is the maximum of ({''') i '''m is the 
time when C, = Cm, and cr is the sign oi ip td, t = tm- 
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FIG. 2: (color online) (a) A typical C, solution (blue solid) 
when A<j!>o is in the fractal region, and the corresponding M[t) 
curve (red dashed); (b) plots of the perturbed solution C,(t) 
(blue solid), the unperturbed solution (red dashed), as well 
as the separatrix solution (black dash-dotted) from the n-th 
saddle approach to the next. 



In the general case where e ^ 0, Eq. ^ is still a 
Hamiltonian system with the conserved Hamiltonian 



i7(C,C,^,^) = i? + £^V[2(l + e)], 



(10) 



where E is given in Eq. ([6]). But E,M and Im(F) are 
not conserved anymore. In order to determine the fractal 
structures as shown in Fig. 1, we need to calculate the 
exit velocity, which corresponds to \C,oo\ for Eq. If 
the orbit has non-zero exit velocity, then from Eqs. ([6]) 
and (IIOD. we find that 



ICool - ^H+^/lPTAPJ{l+e). (11) 

Since H is conserved, to get |Coo|, we only need to find 
Moo- For arbitrary values of e, it is impossible to calcu- 
late Moo analytically. However, when e <C 1 as in Fig. 1 
(where e = 0.001 for both nonlinearities) , the calculation 
of Moo can be done. In this case, Eq. ^ is weakly per- 
turbed from the integrable case (e = 0), thus we will use 
asymptotic techniques in our calculations below. 

To motivate our analysis, we first illustrate in Fig. 2(a) 
a typical C, solution when the initial condition lies in the 
sensitive region of Fig. 1. We see that C undergoes several 
large oscillations, then escapes to — oo. Each oscillation 
corresponds to a "bounce" in the two-wave interactions. 
These bouncing sequences are the key to the existence 
of fractal structures, similar to other physical systems 
IE [1, S B 111 ■ Each local minimum of C will be called a 
saddle approach [ll|. The corresponding M{t) curve is 
also plotted in Fig. 2(a). We see that M changes very 
little near a saddle approach, but changes significantly 
near maxima of C- Below we will calculate the change in 
M from one saddle approach to another. It turns out the 
M formula will be coupled to E and F, thus we need to 
calculate the changes in i?, M and F simultaneously. To 
carry out these calculations, we notice two facts. One is 
that from one saddle approach to the next, the perturbed 
and unperturbed (i.e. integrable) solutions remain close 
to each other since £ <C 1. The other fact is that at 
each saddle approach, £'„,Af„ <C 1. This is so since ini- 
tial i?o, Mq are always small for weak wave interactions, 
and they will remain small when e <C 1. To simplify our 



3 



analysis, we also assume that at each saddle approach, 
M„/En ^ 1. This assumption is satisfied for many ini- 
tial conditions such as 

Now we calculate E, M and F from the n-th saddle 
approach to the next. From Eqs. (|4]) and (O, we get 



AM„ = / —dr 



ee'^ sintpCdT. (12) 



In view of the first fact above, the perturbed orbit in the 
above formula can be approximated by the integrable 
orbit [see Fig. 2(b)]. Due to the second fact, we can 
further approximate the integrable orbit by a separatrix 
orbit. Notice that the separatrix orbit © has three pa- 
rameters. To select the appropriate parameters in the 
separatrix, we note that most contributions to the inte- 
gral of (jl2p come from the ^-maximum region, thus it is 
natural to ask the separatrix solution to have the same C- 
maximum point as the integrable solution [see Fig. 2(b)]. 
Then in view of Eq. ([H]) , the above requirement selects a 
and A in the separatrix ([9]) as 

a„ =. sgn {Im [-2C„coth(C„(T: + F„))]} ; (13) 
A„ = 2 I C^csch^ [C„ « + F„ )] I . (14) 

Here r„ -I- r* is the time the unperturbed solution 
reaches the maximum. Due to the assumption Mn/En <C 
1, the unperturbed solution (to leading order) is a pe- 
riodic solution with period Tp = \pl'K j ^J\E^\. Thus 
T* = rp/2, and r„+i = t„ -|- Tp. Utilizing the above 
results and noticing Tp S> 1, Eq. (fT2|l is asymptotically 
approximated by 



r+oo 



e^' sinipsCsdr = (T„e7ryl„/2. (15) 



By similar calculations and utilizing the symmetry prop- 
erties of the separatrix solution © , we find that 



at;,, 



+ CXD 



(16) 



To calculate Fn+i, notice from Eqs. (H]) and ([7]) that 
F satisfies a linear inhomogeneous ODE: 



dF _ E + iM 
d7 ~ ~ 2(T + iM) 



F + D, 



(17) 



where D is a function of ((^, t}}) whose expression is easy 
to obtain. The homogeneous solution of this ODE is 
C~^{t). Thus by using the method of variation of pa- 
rameters, we can integrate the inhomogeneous ODE (fT7|) 
from T„ to Tn+i and get 



n r"+^ ^E + iMDdT 
Fn+i = T„-f^ + ^^-J ... (18) 



' Cn+l ^JE, 



'ri+1 



- iA/„+i 



Due to the second fact of T„, Af„ <C 1, we can approxi- 
mate the solution (C, "0) in D by the separatrix solution 



Then to leading order in e, we get 



Fn+1 — Fn 



{t + ia)\/E + iM 



n+l 



'n+l 



(19) 



Here t = r — t„ — Tp/2, and a = v^Cn^n^^^- 

Iteration equation is quite complicated. Below, 
we simplify it. From Eq. we get T„ = Tq. Under 

our assumption of A/„/T„ <C 1, to leading order, Eq. 
(Uni) becomes 



Fn+l — 



r)(l- 



M, 



n+l 



At the initial saddle approach, we find from Eq. ([7]) that 



Re(To) = -tt/^2\Eo\. Then solving Eq. we get 

T„ = ^(2|To|)-i/2 [2n-l~ i^„A^„/2To] , (21) 

Sn)Mn- Now we 



where S'„+iA'/„_|_i ~ 2nAf„_|_i — (2n 
introduce a new variable Qn'. 



Qr. 



2nM„ 



-SnMn = -{2\Ea\f'Hm[F,,]/^. (22 



Then substituting formula (PT|) into pS)) . keeping 
only their leading order terms in Mn/En, and putting the 
resulting expressions into (|15p . we obtain the simplified 
iteration equations as 

Mn+i = Mn - sgn(Q„)8|To|'e/^g^, (23) 

Qn+l=Qn + 2Mn+l. (24) 

These equations are derived asymptotically near the sep- 
aratrix orbit ([9]), and will be called the separatrix map. 
This map can be further normalized. Let 

G ^ 8\Eofe/TT, Mn = G-^^Hln, Qn = G'^^^Qn, (25) 

then the normalized separatrix map is 

sgn(e(3„) 



Mn+l = A'ln 
Qn+1 = Qn 



Ql 

2Mn+l. 



(26) 
(27) 



This is a simple but important second-order area- 
preserving map, and it does not have any parameters in 
it (except a sign of e) . This universal map governs weak 
two- wave interactions in generalized NLS equations ([T]). 

Now we compare this map's predictions with direct 
PDE simulations. Here we take the cubic-quintic non- 
linearity in with 7 and initial conditions as for Fig. 
1(a). For the map, we iterate it to infinity to get M^o (in 
practice, 500 iterations performed), which in turn gives 
I Coo I from formula pT|) . With the variable scahngs ([5]) 
and ((25)) considered, the exit velocity graph predicted 
from the map (l^^ - (P7| is shown in Fig. 3(b), while that 
from the PDE simulations is shown in Fig. 3(a). Com- 
paring these two graphs, it is clear that the map gives 
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FIG. 3: (color online) (a,b) The exit-velocity versus initial- 
phase-difference graphs from (a) PDE simulations and (b) 
map predictions; (c) the map's |Moo| graph in the {Qo,Mq) 
plane; (d,e) exit-velocity versus initial-phase-diflterence graphs 
of the PDE for e = 0.01 and 0.001 respectively; (f) plots 
of the accumulation points <jf>L and (I)r versus e from PDE 
simulations (circles and squares) and map predictions (solid 
and dashed lines). 



a good replication of the PDE's fractal structure both 
quahtatively and quantitatively. 

The separatrix map (P5|) - (P7)) exhibits a fractal struc- 
ture in the graph of |Moo| as a function of initial val- 
ues {Qo,Mo), which is displayed in Fig. 3(c) (with 
sgn(e) = 1). This fractal of the map completely de- 
termines the fractal structures in the PDEs For in- 
stance, for initial conditions ([3]), the corresponding initial 
values of the map are 

Qo = Mo = -2-i(7r/e)i/3tan(A0o). (28) 
As A0O varies, Eq. (I^S)) gives a parameterized curve 



in the {Qo,Mo) plane, which is the white straight line 
in Fig. 3 (c). This line cuts cross the fractal of the 
map in Fig. 3(c), and the intersection is precisely the 
fractal structure as observed in Fig. 1 for the PDE 
[see also Fig. 3(a,b)]. The scaling laws for fractals of 
the PDEs can be readily derived from Eq. and 
Fig. 3(c). Let Pi and P2 denote the two accumula- 
tion points of the map's fractal, which are (Qq,Mq) — 
(1.271,1.271) and (0.741,0.741) as marked in Fig. 3(c). 
The corresponding accumulation points (t>ii and 0l in 
the fractal structures of the PDEs are marked in Fig. 
3(d). Here 0l and 0^ are the left and right ends of 
the fractal region. Then according to scalings (|28p . 
we find that (j)R = -ata.n{lA82e^/^ /n^/^), and 0l = 
— atan(2.542e-'^/^/7r^/'^). Hence the map analytically pre- 
dicts that the fractal region of the PDE shrinks to 
A0O = as e 0, and its shrinking speed is propor- 
tional to e^^^ for £ ^ 1. This is precisely what happens. 
To illustrate, we choose two 7 values 0.0029 and 0.0003 
in Eq. which correspond to £ = 0.01 and 0.001 re- 
spectively. The fractal structures of the PDEs for these 
7 values are displayed in Fig. 3(d,e). It is seen that the 
fractal region indeed shrinks as e decreases. We further 
recorded the (f>L and (pR values in the PDE fractals at a 
number of £ values, and the data is plotted in Fig. 3(f). 
The theoretical formulae of (pL and (jjR above are also 
plotted for comparison. It is seen that the PDE values 
and the map's analytical predictions agree perfectly, con- 
firming the asymptotic accuracy of the map (P5|) - (|27p . 

In summary, we have asymptotically analyzed weak 
solitary wave interactions in the generalized nonlinear 
Schrodinger equations and obtained a simple universal 
map. This map gives a complete analytical characteriza- 
tion of universal fractal structures in these wave interac- 
tions. We expect that this work will stimulate research 
in other physical systems where weak solitary wave inter- 
actions arise, such as nonlinear optics, water waves and 
Bose-Einstein condensates. 
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